##This needs MORE MORE WORK
##What to do about sigma0? Is this right.
##Need to add a draw, someting like, fatal[i] ~ dnorm(mu[i],sigma0)

 model{
	for(i in 1:N) { ##observation
		##form the linear predictor (No intercept and constrained Thresholds)
		mu[i] <- 	x[i,1]*beta[1] +  ##SPII
					x[i,2]*beta[2] +  ##caps
					x[i,3]*beta[3] +  ##Majpow
					x[i,4]*beta[4] +  ##Threat
					x[i,5]*beta[5] +  ##Mids
					x[i,6]*beta[6] +  ##Counter
##					x[i,7]*beta[7] +  ##US inv
					alpha1[namehigh[i]] -	 ##Random Effect for Dem
					alpha2[namelow[i]]		 ##Random Effect for Non-Dem
		ystar[i] ~ dnorm(mu[i],prec0)
 		##Cum. Probit Probs.
		##Cat 1
		Q[i,1] <- phi((tau[1] -ystar[i])/sigma0) #tau[1] is the first cutpoint 
		p[i,1] <- Q[i,1]
		
		##Cat 2
		Q[i,2] <- phi((tau[2] -ystar[i])/sigma0) #tau[2] is the second cutpoint
		##slice cdf
		p[i,2] <- Q[i,2]-Q[i,1]
		
		##Cat 3
		Q[i,3] <- phi((tau[3] -ystar[i])/sigma0) #and so on....
		##slice cdf
		p[i,3] <- Q[i,3]-Q[i,2]
		
		##Cat 4
		Q[i,4] <- phi((tau[4] - ystar[i])/sigma0) 
		##slice cdf
		p[i,4] <- Q[i,4]-Q[i,3]
		
		##Cat 5
		Q[i,5] <- phi((tau[5] - ystar[i])/sigma0)
		##slice cdf
		p[i,5] <- Q[i,5]-Q[i,4]
		
		##Cat 6
		Q[i,6] <- phi((tau[6] - ystar[i])/sigma0)
		##slice cdf
		p[i,6] <- Q[i,6]-Q[i,5]
		
		##Cat 7
		Q[i,7] <- phi((tau[7] - ystar[i])/sigma0)
		##slice cdf
		p[i,7] <- Q[i,7]-Q[i,6]
		
		##Cat 8
		Q[i,8] <- phi((tau[8] - ystar[i])/sigma0)
		##slice cdf
		p[i,8] <- Q[i,8]-Q[i,7]
		
		##Cat 9
		Q[i,9] <- phi((tau[9] - ystar[i])/sigma0)
		##slice cdf
		p[i,9] <- Q[i,9]-Q[i,8]
		
		##Cat 10
		Q[i,10] <- phi((tau[10] - ystar[i])/sigma0)
		##slice cdf
		p[i,10] <- Q[i,10]-Q[i,9]
		
		##Cat 11
		p[i,11] <- 1-Q[i,10]
		y[i] ~ dcat(p[i,1:11])## p[i,] sums to 1 for each i
	}
	## Random Effects
	for(p1 in 1:J1) {
	alpha1[p1] ~ dnorm(mu1,prec1)
	}
	for(p2 in 1:J2) {
	alpha2[p2] ~ dnorm(mu2,prec2)
	}
	
	##Prior for sigma0
	prec0<-pow(sigma0,-2)
	sigma0 ~ dunif(0,10000)
	
	##Priors over random effects
	mu1 ~ dnorm(0,.0001)
	mu2 ~ dnorm(0,.0001)
    prec1 <- pow(sigtemp1, -2)
	prec2 <- pow(sigtemp2, -2)
	sigtemp1 ~ dunif(0,10000)
	sigtemp2 ~ dunif(0,10000)

	##priors over beta
	beta[1:6] ~ dmnorm(b0[],B0[,])
	
}


